The paper Sampling and homology via bottlenecks by Di Rocco et. al. introduces a novel algorithm for producing a provably dense sampling of a smooth and compact algebraic variety where the density is determined by the smallest bottleneck. Using the sample, the zeroth and first homology of the variety can be recovered using a modified Vietoris-Rips construction. This notebook implements the algorithm for sampling and homology computation for the case of complete intersections and illustrates it using a curve in $\Bbb R^2$.
using HomotopyContinuation, DynamicPolynomials, LinearAlgebra, IterTools, Random
Random.seed!(1)
n=5 # ambient dimension
@polyvar x[1:n] y[1:n] p[1:n] gamma[1:n]
singular_cubic = x[1]*(x[2]^2+x[3]^2+x[4]^2+x[5]^2)-(x[2]^3+x[3]^3+x[4]^3-(1/2)*x[5]^3);
smoothing = singular_cubic-(1/4)*x[1]^3;
deg_three_generic = rand((-1,1)) + rand((-1,1))*x[1] + rand((-1,1))*x[1]^2 + rand((-1,1))*x[1]^3 + rand((-1,1))*x[2] + rand((-1,1))*x[1]*x[2] + rand((-1,1))*x[1]^2*x[2] + rand((-1,1))*x[2]^2 + rand((-1,1))*x[1]*x[2]^2 + rand((-1,1))*x[2]^3 + rand((-1,1))*x[3] + rand((-1,1))*x[1]*x[3] + rand((-1,1))*x[1]^2*x[3] + rand((-1,1))*x[2]*x[3] + rand((-1,1))*x[1]*x[2]*x[3] + rand((-1,1))*x[2]^2*x[3] + rand((-1,1))*x[3]^2 + rand((-1,1))*x[1]*x[3]^2 + rand((-1,1))*x[2]*x[3]^2 + rand((-1,1))*x[3]^3 + rand((-1,1))*x[4] + rand((-1,1))*x[1]*x[4] + rand((-1,1))*x[1]^2*x[4] + rand((-1,1))*x[2]*x[4] + rand((-1,1))*x[1]*x[2]*x[4] + rand((-1,1))*x[2]^2*x[4] + rand((-1,1))*x[3]*x[4] + rand((-1,1))*x[1]*x[3]*x[4] + rand((-1,1))*x[2]*x[3]*x[4] + rand((-1,1))*x[3]^2*x[4] + rand((-1,1))*x[4]^2 + rand((-1,1))*x[1]*x[4]^2 + rand((-1,1))*x[2]*x[4]^2 + rand((-1,1))*x[3]*x[4]^2 + rand((-1,1))*x[4]^3 + rand((-1,1))*x[5] + rand((-1,1))*x[1]*x[5] + rand((-1,1))*x[1]^2*x[5] + rand((-1,1))*x[2]*x[5] + rand((-1,1))*x[1]*x[2]*x[5] + rand((-1,1))*x[2]^2*x[5] + rand((-1,1))*x[3]*x[5] + rand((-1,1))*x[1]*x[3]*x[5] + rand((-1,1))*x[2]*x[3]*x[5] + rand((-1,1))*x[3]^2*x[5] + rand((-1,1))*x[4]*x[5] + rand((-1,1))*x[1]*x[4]*x[5] + rand((-1,1))*x[2]*x[4]*x[5] + rand((-1,1))*x[3]*x[4]*x[5] + rand((-1,1))*x[4]^2*x[5] + rand((-1,1))*x[5]^2 + rand((-1,1))*x[1]*x[5]^2 + rand((-1,1))*x[2]*x[5]^2 + rand((-1,1))*x[3]*x[5]^2 + rand((-1,1))*x[4]*x[5]^2 + rand((-1,1))*x[5]^3;
deg_two_generic = rand((-1,1)) + rand((-1,1))*x[1] + rand((-1,1))*x[1]^2 + rand((-1,1))*x[2] + rand((-1,1))*x[1]*x[2] + rand((-1,1))*x[2]^2 + rand((-1,1))*x[3] + rand((-1,1))*x[1]*x[3] + rand((-1,1))*x[2]*x[3] + rand((-1,1))*x[3]^2 + rand((-1,1))*x[4] + rand((-1,1))*x[1]*x[4] + rand((-1,1))*x[2]*x[4] + rand((-1,1))*x[3]*x[4] + rand((-1,1))*x[4]^2 + rand((-1,1))*x[5] + rand((-1,1))*x[1]*x[5] + rand((-1,1))*x[2]*x[5] + rand((-1,1))*x[3]*x[5] + rand((-1,1))*x[4]*x[5] + rand((-1,1))*x[5]^2;
epsilon = 1/100
F = [
smoothing+epsilon*deg_three_generic,
x[1]^2+x[2]^2+x[3]^2+x[4]^2+x[5]^2-1+epsilon*deg_two_generic
]
d=length(F) # codimension of variety
k = n-d # dimension of variety
@polyvar lambda[1:d] mu[1:d]; # lagrange multipliers
# Compute the bottlenecks
grad = differentiate(F, x)
G = subs(F, x => y)
grady = subs(grad, x => y)
system = [F; G; map(j -> x[j]-y[j]-dot(lambda, grad[:, j]), 1:n); map(j -> x[j]-y[j]-dot(mu, grady[:, j]), 1:n)]
result = solve(system, start_system = :total_degree)
Tracking 2125764 paths... 100%|████████████████████████████████████████| Time: 1:53:40 # paths tracked: 2125764 # non-singular solutions (real): 20724 (712) # singular solutions (real): 5721 (12) # total solutions (real): 26445 (724)
Result{Vector{ComplexF64}} with 23856 solutions
===============================================
• 20724 non-singular solutions (712 real)
• 3132 singular solutions (12 real)
• 2125764 paths tracked
• random seed: 835247
• multiplicity table of singular solutions:
┌───────┬───────┬────────┬────────────┐
│ mult. │ total │ # real │ # non-real │
├───────┼───────┼────────┼────────────┤
│ 1 │ 543 │ 12 │ 531 │
│ 2 │ 2589 │ 0 │ 2589 │
└───────┴───────┴────────┴────────────┘
# Choose the size of the grid as the smallest bottleneck
bottlenecks = map(s -> (s[1:n], s[n+1:2*n]), real_solutions(nonsingular(result)))
bn_lengths = sort!(map(b -> (norm(b[1]-b[2]), b), bottlenecks), by = a -> a[1])
δ = bn_lengths[1][1]/2 # grid size
0.09505817482972509
# Compute the bounding box by computing the EDD starting from the center of the largest bottleneck
q = bn_lengths[end][2][1] + (bn_lengths[end][2][2]-bn_lengths[end][2][1])/2 + randn(Float64, n)*10^(-10)
system = [F; map(j -> x[j]-q[j]-dot(lambda, grad[:, j]), 1:n)]
result = solve(system, start_system = :polyhedral)
Result{Vector{ComplexF64}} with 156 solutions
=============================================
• 156 non-singular solutions (24 real)
• 0 singular solutions (0 real)
• 156 paths tracked
• random seed: 104798
# Extract farthest point from q to X and use as box length
critical_points = sort!(map(c -> (norm(c[1:n]-q), c[1:n]), real_solutions(nonsingular(result))), by = a -> a[1])
b = critical_points[end][1]
indices = [i for i in -b:δ:b];
# Compute basic sample
samples = []
counter = 0
start_time = time_ns()
for s in IterTools.subsets(1:n, k)
Ft = [F; map(i -> x[s[i]]-p[i]-q[s[i]], 1:k)]
p₀ = randn(ComplexF64, k)
F_p₀ = subs(Ft, p[1:k] => p₀)
result_p₀ = solve(F_p₀)
S_p₀ = solutions(result_p₀)
# Construct the PathTracker
tracker = HomotopyContinuation.pathtracker(Ft; parameters=p[1:k], generic_parameters=p₀)
for p1 in Iterators.product(map(j-> 1:length(indices), s)...)
counter += length(S_p₀)
for s1 in S_p₀
result = track(tracker, s1; target_parameters=map(j -> indices[p1[j]], 1:k))
# check that the tracking was successfull
if is_success(result) && is_real(result)
push!(samples, real(solution(result)))
end
end
end
end
time_basic_sample = time_ns()-start_time;
# Compute extra sample
extra_samples = []
extra_counter = 0
start_time = time_ns()
for l in 1:k-1
for s in IterTools.subsets(1:n, l)
Ft = [F; map(i -> x[s[i]]-p[i]-q[s[i]], 1:l)]
grad = differentiate(Ft, x)
system = [Ft; map(j -> x[j]-y[j]-dot(gamma[1:n-k+l], grad[:, j]), 1:n)]
p₀ = randn(ComplexF64, n+l)
F_p₀ = subs(system, [y; p[1:l]] => p₀)
result_p₀ = solve(F_p₀)
S_p₀ = solutions(result_p₀)
# Construct the PathTracker
tracker = HomotopyContinuation.pathtracker(system; parameters=[y; p[1:l]], generic_parameters=p₀)
for p1 in Iterators.product(map(j-> 1:length(indices), s)...)
extra_counter += length(S_p₀)
for s1 in S_p₀
result = track(tracker, s1; target_parameters=[randn(Float64, n); map(j -> indices[p1[j]], 1:l)])
# check that the tracking was successfull
if is_success(result) && is_real(result)
push!(extra_samples, real(solution(result))[1:n])
end
end
end
end
end
time_extra_sample = time_ns()-start_time;
Tracking 1458 paths... 100%|████████████████████████████████████████| Time: 0:00:00 # paths tracked: 1458 # non-singular solutions (real): 24 (0) # singular solutions (real): 0 (0) # total solutions (real): 24 (0)
# Summary
println("grid size: ", δ)
println("length of bounding box: ", b)
println("number of tracked paths: ", counter+extra_counter)
println("total time: ", (time_basic_sample+time_extra_sample)*10^(-9), " s")
println("|E_δ|: " , length(samples), ", |E'_δ|: ", length(extra_samples))
println("total number of samples: ", length(samples)+length(extra_samples))
grid size: 0.09505817482972509 length of bounding box: 1.0137680752023333 number of tracked paths: 762300 total time: 83.85367645500004 s |E_δ|: 69154, |E'_δ|: 10841 total number of samples: 79995
using DelimitedFiles
writedlm( "pointcloud.csv", samples, ',')
using Plots
plotlyjs()
ArgumentError: Package PlotlyJS not found in current path:
- Run `import Pkg; Pkg.add("PlotlyJS")` to install the PlotlyJS package.
Stacktrace:
[1] require(into::Module, mod::Symbol)
@ Base ./loading.jl:967
[2] top-level scope
@ ~/.julia/packages/Plots/hyS17/src/backends.jl:491
[3] eval
@ ./boot.jl:373 [inlined]
[4] _initialize_backend(pkg::Plots.PlotlyJSBackend)
@ Plots ~/.julia/packages/Plots/hyS17/src/backends.jl:490
[5] backend(pkg::Plots.PlotlyJSBackend)
@ Plots ~/.julia/packages/Plots/hyS17/src/backends.jl:174
[6] #plotlyjs#242
@ ~/.julia/packages/Plots/hyS17/src/backends.jl:31 [inlined]
[7] plotlyjs()
@ Plots ~/.julia/packages/Plots/hyS17/src/backends.jl:31
[8] top-level scope
@ In[26]:2
[9] eval
@ ./boot.jl:373 [inlined]
[10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1196
# Plot the sampled points
S = reduce(hcat, vcat(samples, extra_samples))
if n < 3
scatter(S[1,:], S[2,:])
else
scatter(S[1,:], S[2,:], S[3,:])
end
using Eirene
#Compute distance matrix
ϵ = 2*δ
D = zeros((length(samples), length(samples))) #Initialize distance matrix
neighbour_lists = []
for i in 1:length(samples)
push!(neighbour_lists, [])
end
candidate_lists = []
for i in 1:length(samples)
candidate_list = []
for j in (i+1):length(samples)
dist = norm(samples[i]-samples[j])
if dist < sqrt(8)*δ
if dist < ϵ
push!(neighbour_lists[i], j)
push!(neighbour_lists[j], i)
else
push!(candidate_list, j)
end
end
D[i, j] = dist
D[j, i] = dist
end
push!(candidate_lists, candidate_list)
end
OutOfMemoryError()
Stacktrace:
[1] Array
@ ./boot.jl:459 [inlined]
[2] Array
@ ./boot.jl:467 [inlined]
[3] zeros
@ ./array.jl:525 [inlined]
[4] zeros(dims::Tuple{Int64, Int64})
@ Base ./array.jl:522
[5] top-level scope
@ In[29]:4
[6] eval
@ ./boot.jl:373 [inlined]
[7] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1196
# Modify distance matrix to add edges in VR-complex
DM = deepcopy(D) #Modified distance matrix
thresh = 2*δ - 10^(-10)
for i in 1:length(samples)
for j in candidate_lists[i]
for k in neighbour_lists[j]
if D[i, k] < ϵ && D[j, k] < ϵ
DM[i, j] = thresh
DM[j, i] = thresh
break
end
end
end
end
UndefVarError: D not defined Stacktrace: [1] top-level scope @ In[30]:3 [2] eval @ ./boot.jl:373 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1196
# Homology computation using Eirene using VR-complex with distance ϵ
C = eirene(DM, model="vr", maxdim=1, minrad=ϵ, maxrad=ϵ);
UndefVarError: DM not defined Stacktrace: [1] top-level scope @ In[31]:2 [2] eval @ ./boot.jl:373 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1196
println("0-th Betti number: ", betticurve(C, dim=0)[1, 2])
println("1-th Betti number: ", betticurve(C, dim=1)[1, 2])
UndefVarError: C not defined Stacktrace: [1] top-level scope @ In[32]:1 [2] eval @ ./boot.jl:373 [inlined] [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1196